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INTRODUCTION 


The  propagation  of  elastic-plastic  waves  in  long  rods  has  been  treated 
extensively  in  the  literature  (refs  1-3)  since  the  pioneering  works  of  Donell, 
Karraan,  Taylor,  and  Rakhmatulin.  The  study  of  plastic  wave  propagation  is 
important  because  it  attempts  to  explain  the  response  of  materials  to  intense 
dynamic  loading  and  serves  also  as  a  basis  for  determining  dynamic  material 
properties. 

Analytical  solutions  can  be  obtained  for  only  a  few  idealized  situations; 
hence,  many  impact  studies  have  been  performed  using  numerical  methods.  Many 
computer  codes  using  either  finite-element  or  finite-difference  approaches 
have  been  developed.  The  computer  simulation  of  impact  phenomena  in  solids  is 
still  quite  involved  and  it  depends  critically  on  the  impact  velocity.  For 
high  velocity  impact  and  penetration  problems,  a  good  review  was  given  by 
Zukas  et  al  (ref  3).  For  low  velocity  contact-impact  problems,  many 
structural  response  codes  were  reviewed  by  Noor  (ref  4). 

A  numerical  study  of  the  dynamic  response  of  an  elastic-plastic 
projectile  due  to  normal  impact  is  reported  here  using  the  finite  element 
structural  response  code  ADINA  (ref  5).  The  projectile  is  a  finite  length 

^•Cristescu,  N.,  Dynamic  Plasticity,  John  Wiley  &  Sons,  Inc.,  New  York,  1967. 
^Nowacki,  W.  K. ,  Stress  Waves  in  Non-Elastic  Solids,  PerRaraon  Press,  Oxford. 
1978. 

o 

JZukas,  J.  A.,  Nicholas,  T. ,  Swift,  H.  F.,  Greszczak,  L.  B.,  and  Curran,  D. 

R.,  Impact  Dynamics,  John  Wiley  &  Sons,  Inc.,  New  York,  1982. 

^Noor,  A.  K. ,  "Survey  of  Computer  Programs  for  Solution  of  Nonlinear 
Structural  and  Solid  Mechanics  Problems,"  Computer  and  Structures,  Vol.  13, 
1981,  pp.  425-465. 

5Bathe,  K.  J.,  "ADINA  Users'  Manual,"  Report  AE-81-1,  ADINA  Engineering,  Inc., 
Watertown,  MA,  1981. 
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cylindrical  bar  made  of  a  high  strength  steel.  The  bar  is  long  and  travels 
with  velocity  V  =  75  m/s  before  it  strikes  a  rigid  target.  First,  three 
direct  integration  schemes  have  been  used  for  the  uniaxial  stress  wave  problem 
in  a  linear-hardening  material,  and  the  results  are  compared  with  an  exact 
analytical  solution  in  order  to  evaluate  the  accuracy  and  stability.  Then, 
additional  numerical  results  for  perf ectly-plastic  materials  are  discussed  in 
order  to  show  the  effect  of  strain-hardening  for  a  multi-linear  material 
model.  Finally,  some  results  based  on  two-dimensional  elements  are  presented 
in  order  to  show  the  lateral  effect. 

ANALYTICAL  SOLUTION 

The  problem  of  the  normal  impact  of  a  rod  against  a  rigid  target  has  been 
considered  by  many  authors.  Various  schemes  have  been  used  for  different 
kinds  of  initial  conditions  and  various  material  properties.  For  a  linear 
work-hardening  material  due  to  sudden  impact,  an  analytical  solution  for  the 
uniaxial  stress  wave  problem  is  available  (ref  1)  and  it  is  presented  here  for 
comparison  with  the  corresponding  ADINA  results.  Thus,  consider  a  bar  of 
length  L  and  diameter  D  that  is  moving  with  velocity  V  in  the  negative  Z 
direction.  At  time  0  the  bar  strikes  a  rigid  wall  Z  =  0  (Figure  la) .  Guided 
by  the  stress-strain  curve  for  a  high  strength  steel  supplied  to  us  (ref  6), 
the  following  material  data  will  be  used: 

E  =  208  GPa,  p  =  0.783  g/cc,  v  =  0.293 

(1) 

=  1.3  GPa,  Ep  =  4  GPa, 

^Cristescu,  N.  ,  Dynamic  Plasticity,  John  Wiley  &  Sons,  Inc.,  New  York,  1967. 
^Wright  T.  W . ,  Private  Communication,  1982. 
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where  E,  p,  v,  oy,  and  Ep  are  Young’s  modulus,  density,  Poisson’s  ratio,  yield 
stress,  and  plastic  modulus,  respectively.  The  elastic  and  plastic  one¬ 


dimensional  wave  speeds  will  then  be  ce  =  vE/p  =  5154  m/s  and  cp  =  /Ep/ p  - 
715  m/s,  respectively.  The  velocity  Vy  corresponding  to  the  yield  stress  oy 
of  the  material  is 

Vy  =  oy/(pce)  -  ceCy  =  32.21  m/s  (2) 

Both  elastic  and  plastic  waves  will  be  generated  if  the  impact  velocity  V  > 

Vy.  Let  us  consider  the  case  V  =  75  m/s,  L  =  1.1  m,  D  =  0.1  m.  After  impact, 
two  shock  wave  fronts  delimit  three  distinct  regions  in  the  bar  (Figures  lb 
and  lc) .  The  analytical  solutions  for  the  particle  velocity,  strain,  and 
stress  in  these  three  regions  are 

V o  =  “V  =  -75  m/s,  =0,  0O  =  0, 

Vf  =  -V  +  Vy  =  -42.79  ra/s,  =  -£y  -  -0.625%,  aQ  =  -oy  =  -1.3  GPa, 

V2  =  0,  e2  =  -ty  +  (vy-v)/cp  =  -6.610% 

°2  =  -^y  +  Ep(£2+£y)  =  -1.539  GPa  .  (3) 

After  time  t  =  L/ce  *  213  us,  the  elastic  wave  front  is  reflected  from  the 
free  end.  Behind  this  front  (Figure  Id),  03  =  e3  =  0,  V3  =  2Vy  -  V.  At  time 
tg  =  2L/(ce+cp)  =  374.85  us,  the  wave  fronts  of  the  plastic  and  of  the  elastic 
unloading  waves  meet  at  the  section  S.  The  stress  and  velocity  are  continuous 
but  the  strain  is  discontinuous  across  this  section  S,  a  nonpropagable 
discontinuity  surface.  Since  the  inequality 

2ce 

2Vv  <  V  <  (1  +  - )  Vv  =  88.79  m/s  (4) 

ce+cp 

is  satisfied,  the  plastic  wave  stops  at  S  and  elastic  waves  again  propagate 
from  S  in  both  directions  (Figure  le) .  The  analytical  solutions  in  regions  4 
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and  5  are 


1  ce 

V4  =  V5  =  -  (3  - )(Vy-V)  +  V  =  13.79  ra/s 

2  Cp 

1 

G4  =  -  (l+Cp/ce) ( £y-V/ce)  =  -0.473% 
c  c 

£5  =  -  (2  --  +  1  -  — — ) ( Ey— V / c e )  =  -6.342% 
cp  ce 

a4  =  a5  =  -0.983  GPa  (5) 

Figures  lb  through  le  show  the  locations  of  the  wave  fronts  at  time  t  = 
100,  200,  300,  400  ps,  respectively.  The  analysis  can  be  continued  until  the 
contact  between  the  bar  and  the  target  ceases  at  tc  =  4L/(ce+Cp)  =  749.7  ps . 
More  detailed  information  about  the  analytical  solution  can  be  found  in  the 
book  by  Gristescu  (ref  1). 

ADINA  SOLUTION 

The  ADINA  code,  developed  by  K.  J.  Bathe,  is  a  general  purpose  finite 
element  program  for  Automatic  Dynamic  Incremental  Nonlinear  Analysis  (ref  5). 
In  nonlinear  analysis  the  incremental  finite  element  equations  of  motion  used 
are,  in  implicit  time  integration, 

M  t+Aty  +  c  t+Atu  +  CK  U  =  t+AtR  -  c£  (6) 

and  in  explicit  time  integration, 

M  CU  +  _C  CU  +  c K  U  =  CR  -  C_F  (7 ) 

where  M,  C,  CK,  CR,  t+AtR,  CF  are  constant  mass  matrix,  constant  damping 
matrix,  tangent  stiffness  matrix  at  time  t,  external  load  vector  applied  at 

^Gristescu,  N. ,  Dynamic  Plasticity,  John  Wiley  &  Sons,  Inc.,  New  York,  1967. 
^Bathe,  K.  J. ,  "ADINA  Users1  Manual,"  Report  AE-81-1,  ADINA  Engineering,  Inc., 
Watertown,  MA,  1981. 
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time  t,  t+At,  nodal  point  force  vector  equivalent  to  the  element  stresses  at 


time  t,  respectively,  and  IJ  is  the  vector  of  nodal  point  displacement 
increments  from  time  t  to  time  t+At,  i.e.,  IJ  =  t+Aty  -  tjj,  The  solution  of 
Eq.  (6)  yields,  in  general,  an  approximate  displacement  increment  JJ.  To 
improve  the  solution  accuracy  and  in  some  cases  to  prevent  the  development  of 
numerical  instabilities,  it  may  be  necessary  to  use  equilibrium  iteration  in 
each  or  preselected  time  steps. 

In  ADINA,  the  central  difference  method  is  employed  for  explicit  time 
integration  and  either  the  Newmark  method  or  Wilson  method  is  employed  for 
implicit  time  integration.  The  integration  schemes  (ref  7)  are  given  by: 

CU  =  -----  (t+Atjj  _  2tu  +  t-Atu)  (8 

(At)  2 

•  1 

tu  - - (t+Atu  -  t-Atu)  (9 

2  At 

for  the  central  difference  method, 


(ii) 


(10) 


for  the  Newmark  method,  and 


t+Atu  =  tu  +  (---) (t+0Atu  -  tu) 


(12) 


e  >  1,  0  <  t  <  0At 


for  the  Wilson  method. 


7Bathe,  K.  J.,  Numerical  Methods  in  Finite  Element  Analysis,  Prentice-Hall, 
Inc.,  New  Jersey,  1976. 


5 


The  Wilson  and  Newmark  methods  are  unconditionally  stable  if  0  >  1.37  or 
a  >  l/4(l/2+6),  6  >  1/2.  In  our  numerical  study,  we  have  chosen  0  =  1*4,  a  = 
1/4,  6  =  1/2.  In  using  the  central  difference  method,  the  time  step,  At,  has 
to  satisfy  the  Courant  condition 

2 

At  =  Atcr  =  KA£/C  or  -  (13) 

0) 

where  A£  is  the  minimum  mesh  size,  u)  is  the  maximum  natural  frequency,  c  is 
the  local  sound  speed,  and  K  <  1. 

NUMERICAL  COMPARISON 

Consider  a  bar  with  the  following  geometrical  and  material  data:  L  - 
1.1  m,  D  -  0.1  m,  E  =  208  GPa,  p  =  0.783  g/cc,  ay  =  1.3  GPa,  Ep  =  4  GPa, 
subjected  to  two  values  of  impact  velocity:  V  =25  m/s  or  75  m/s.  Since  the 
velocity  corresponding  to  the  yield  stress  Oy  of  the  material  is  Vy  =  32.2 
m/s,  the  impact  is  elastic  for  the  first  case  and  elastic-plastic  for  the 
second  case.  Analytical  solutions  are  known  for  both  cases.  We  used  100 
one-dimensional  truss  elements  to  simulate  this  uniaxial  stress  wave  problem. 
In  order  to  satisfy  the  stability  criterion  for  explicit  integration  by  the 
central  difference  method,  we  have  chosen  the  time  step  At  =  2  ps  which  is  less 
than  the  critical  time  step 

Atcr  =  AZ/ce  =  2.13  ps 

Our  study  compares  three  integration  schemes  with  the  same  time  step. 

The  computations  are  all  stable,  and  the  numerical  results  for  the  axial 
stress  and  velocity  when  V  =  25  m/s  are  shown  in  Figures  2  through  5.  Figure 
2  shows  results  for  the  particle  velocity  and  stress  along  the  rod  at  t  =  100, 
200,  300,  400,  500  ps  when  the  Wilson  method  was  used.  Figure  3  shows  the 
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similar  results  for  the  Newmark  method.  The  numerical  results  based  on  these 


two  methods  are  less  accurate  when  compared  with  the  results  based  on  the 
central  difference  method.  Figure  4  shows  a  comparison  of  the  results  for  the 
particle  velocity  along  the  length  of  the  bar  based  on  the  central  difference 
method,  the  Newmark  method,  and  the  analytical  solution  at  t  =  100,  200,  300, 
400  ps.  A  similar  comparison  between  the  central  difference  method,  and  the 
analytic  solution  for  the  axial  stress  along  the  bar  is  shown  in  Figure  5.  It 
should  be  noted  that  the  computations  were  performed  under  the  assumption  that 
the  rod  and  target  remain  in  contact  after  impact  while  the  theoretical 
interval  of  contact  is  tc  =  2L/ce  =  426  ps. 

Calculations  were  also  performed  for  elastic-plastic  impact  with  V  =  75 
m/s  using  the  three  integration  schemes  with  the  same  mesh  size  and  time  step. 
A  comparison  of  numerical  results  for  the  axial  stress  and  velocity  using  the 
central  difference  and  Newmark  methods  is  shown  in  Figure  6  at  t  =  200  ps . 

The  solid  and  dotted  curves  represent  the  results  based  on  the  central 
difference  method  and  Newmark  method,  respectively.  Similar  results  are  shown 
in  Figure  7  at  t  5  400  ps.  As  can  be  seen  from  these  two  figures,  the  central 
difference  method  gives  more  accurate  results  than  the  Newmark  method.  The 
numerical  results  based  on  the  Wilson  method  are  compared  with  those  based  on 
the  central  difference  method  in  Figure  8.  This  also  demonstrates  that  the 
numerical  results  based  on  the  central  difference  method  are  more  accurate. 

We  may  then  conclude  from  this  study  on  elastic  as  well  as  elastic-plastic 
impact  that  the  central  difference  method  will  give  more  accurate  results  than 
the  other  two  integration  schemes. 
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HARDENING  AND  LATERAL  EFFECTS 


After  reaching  the  above  conclusion,  we  used  the  central  difference 
method  for  the  rest  of  this  study.  In  order  to  show  hardening  effects,  we 
obtained  numerical  results  for  displacement,  velocity,  strain,  and  stress  in  a 
long  rod  of  an  elastic-perf ectly-plastic  material.  Figures  9  and  10  show  the 
results  of  the  particle  velocity  and  stress  for  a  linear  work-hardening  as 
well  as  a  perf ectly-plastic  material  at  t  =  300  and  400  ps,  respectively.  It 
can  be  seen  from  these  comparsions  that  the  effect  of  strain-hardening  on  the 
particle  velocity  and  stress  is  quite  significant  even  though  the  plastic 
modulus  Ep  =  4  GPa  is  small  when  compared  with  the  elastic  modulus  E  =  208 
GPa. 

In  order  to  study  lateral  effects,  we  have  used  two-dimensional  four-node 
quadrilateral  ring  elements  to  obtain  numerical  results.  We  chose  the  same 
mesh  size  Ar  -  Az  =  0.011  m  and  used  50  elements  along  the  length  of  the  bar 
with  L/D  =25.  In  this  arrangement,  the  new  length  of  the  bar  is  only  half 
of  the  original.  We  have  used  the  same  time  step  At  =  2  ps  as  in  the  one¬ 
dimensional  truss  elements.  This  time  step  yields  stable  computations  for 
one-dimensional  truss  elements  but  not  for  two-dimensional  quadrilateral  ring 
elements.  This  seems  due  to  the  lateral  effect  such  as  the  Poisson’s  ratio- 
including  the  effect  of  Poisson’s  ratio  ( v  =  0.293),  the  speed  of  longitudinal 
elastic  wave  is  c<j  =  [ (E/p) ( (l-v)/(l+v))/(l-2  v)  ] =  5923  m/s,  which  reduces 
to  ce  =  (E/p)!/^  =.  5154  m/s  in  case  of  v  =  0.  Guided  by  the  stability 
criterion  for  linear  elastic  problems,  we  thus  chose  At  <  Atcr  =  AZ/c^  = 
0.011/5923  sec  =  1.857  ps.  For  this  reason,  we  carried  out  the  computations 
for  250  time  steps  using  At  =  1  ps.  The  total  number  of  time  steps  is  the 
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same  for  both  the  one  and  two-dimensional  problems.  The  length  of  the  bar  and 
time  increment  for  the  two-dimensional  case  are  only  half  of  the  one¬ 
dimensional  case.  For  the  two-dimensional  case,  we  have  carried  out  the 
computations  for  impact  velocities  of  V  =  25  m/s  and  75  m/s.  The  results  for 
the  elastic  impact  (V  =  25  m/s)  are  not  shown  here.  For  elastic-plastic 
impact  (V  =  75  m/s),  we  have  carried  out  the  computations  for  a  bilinear  as 
well  as  a  multi-linear  material  model.  Seven  points  are  used  to  represent  the 
stress-strain  curve,  i.e.,  (a  in  GPa,  e  in  %)  =  (1.3,  0.63),  (1.355,  1.03), 
(1.38,  1.83),  (1.394,  2.63),  (1.415,  4.23),  (1.43,  5.83),  (1.45,  8.63). 

The  numerical  results  for  the  case  of  a  multi-linear  material  are  shown 
in  Figures  11  through  14.  Figure  11  shows  the  effective  stress  and  axial 
velocity  along  the  length  of  the  bar  at  the  end  of  50  and  100  time  steps. 
Curves  .1  and  2  represent  the  results  at  t  -  50  and  100  ps,  respectively. 

Figure  12  shows  the  axial  stresses  along  the  length  of  the  bar.  We  have  2x2 
stations  to  carry  out  the  numerical  integration  in  the  spatial  direction  in 
each  element.  The  results  for  the  stresses  are  calculated  at  four  integration 
stations.  As  can  be  expected  from  a  long  rod  (L/D  =  25),  the  differences  in 
the  results  along  the  stations  near  the  centerline  or  outside  are  very  small. 
Figures  13  and  14  show  the  similar  results  for  the  effective  stress,  particle 
velocity,  and  axial  stress  along  the  length  of  the  bar  at  t  =  150  and  200  ps. 
As  can  be  seen  in  these  two  figures,  the  axial  stress  in  the  plastic  zone 
shows  bigger  oscillations  than  corresponding  effective  stress.  Comparing  the 
results  shown  in  Figures  11  through  14  with  the  one-dimensional  results  shown 
in  Figures  6  through  8,  we  conclude  that  the  transition  near  the  wave  front  is 
not  as  steep  as  the  one-dimensional  case  and  dispersion  behind  the  wave  front 
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can  be  observed.  Some  of  the  oscillations  are  real  due  to  lateral  effects 
such  as  radial  inertia,  radial  shear,  etc.,  but  some  are  due  to  numerical 
errors  such  as  truncation  error  in  the  finite  element  system,  approximations 
in  time  integration  schemes,  numerical  integration  in  the  spatial  directions, 
etc.  It  is  difficult  to  identify  how  many  of  the  oscillations  are  real  and 
how  many  are  due  to  numerical  error.  We  hope  to  develop  a  numerical  model  to 
minimize  the  numerical  error  for  this  purpose  while  trying  to  improve  the 
theoretical  model. 

CONCLUSION 

Based  on  our  numerical  study  of  the  uniaxial  stress  wave  problem  in  a 
linear-hardening  material,  the  central  difference  method  gives  more  accurate 
results  than  the  Wilson  and  Newmark  methods.  The  effect  of  hardening  is 
significant  and  lateral  effects  due  to  radial  motion  need  further  study.  We 
plan  to  improve  the  theoretical  model  and  to  develop  a  numerical  scheme.  We 
hope  to  compare  our  numerical  results  with  experiments  involving  normal  impact 


of  cylindrical  rods. 
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Figure  1.  Analytical  solution  to  a  projectile  striking  a  rigid  target 
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Figure  2.  Axial  velocity  and  stress  based  on  the  Wilson  method  (V  =  25  m/s). 
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Figure  3.  Axial  velocity  and  stress  based  on  the  Nevmark  method  (V  =  25  m/s). 
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Figure  4.  A  comparison  of  the  central  difference  method,  the  Newmark 
method,  and  the  analytical  solution  for  the  particle 
velocity  (V  =  25  m/s). 
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Figure  5.  A  comparison  of  the  central  difference  method  and  the 
analytical  solution  for  the  axial  stress  (V  =  25  m/s) 
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Figure  6. 


A  comparison  of  the  central  difference  method  and  the  Newmark 
method  for  V  and  a  at  t  =  200  ps  (V  =  75  m/s). 
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Figure  7.  A  comparison  of  the  central  difference  method  and  the  Newmark 
method  for  V  and  a  at  t  =  400  ps  (V  =  75  m/s). 
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Figure  8.  A  comparison  of  the  central  difference  method  and  the  Wilson 
method  for  V  at  t  =  200  and  400  ys  (V  =  75  m/s) . 
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Figure  11.  The  effective  stress  and  axial  velocity  based  on  2-D 
elements  at  t  =  50  and  100  ys. 
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Figure  12. 


Axial  stress  (az)  based  on  2-D  elements  at  t  =  50  and  100  ys. 
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Figure  13.  The  effective  stress  and  axial  velocity  based  on  2-D 
elements  at  t  =  150  and  200  ys. 
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Figure  14.  Axial  stress  (a2)  based  on  2-D  elements  at  t  =  150  and  200  pis. 
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